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THE NASA AMES RESEARCH CENTER ONE- AND TWO-DIMENSIONAL STRATOSPHERIC MODELS 

PART I: THE ONE-DIMENSIONAL MODEL 

R. P. Turco 
R and D Associates 
and 

R. C, Whitten 
Ames Research Center 

SUMMARY 


A one~dimensional model of stratospheric trace constituents, developed 
in a joint effort by scientists at Ames Research Center and at R and D 
Associates, is described in detail. Specifically, the numerical solution 
of the species continuity equations, including a technique for treating 
the "stiff" differential equations representing the chemical kinetic terms, 
and an appropriate method for simulating the diurnal variation of the species 
concentrations, are discussed. A specialized treatment of atmospheric 
photodissociation rates is outlined in the text. The choice of a vertical 
eddy diffusivity profile and its success in predicting the vertical tracer 
distributions (carbon 14, methane, and nitrous oxide) are also discussed. 


INTRODUCTION 


In this report we describe a one-dimensional atmospheric computer 
model. Although many of the technical details are presented here, a deeper 
insight into the model can be obtained from several specialized papers: 
the numerical integration scheme (ref. 1), the calculation of photodissocia- 
tion rates (ref. 2), and the treatment of diurnal averaging (ref. 3). To be 
concise in this report, we emphasize the structure and philosophy of the 
model, and we avoid detailed tables of physical data such as photodissociation 
cross sections which do not affect the model performance. We also stress 
those characteristics of the model that differentiate it from other models. 

To complement this review we refer the reader to several publications that 
illustrate specific model simulations of: carbon compounds in the strato- 

sphere and mesosphere (ref. 4); supersonic transport exhaust contamination 
at high altitudes (ref. 5); the diurnal variations of hydrogen and nitrogen 
compounds in the stratosphere (ref. 6); the environmental impact of hydrogen 
chloride released during space shuttle flights (ref. 7) ; the effects on ozone 
of nitrogen oxides generated by a nuclear war (ref. 8) ; fluorocarbon deple- 
tions of global ozone (ref. 9), and the formation and evolution of strato- 
spheric aerosol particles (ref. 10). 



SPECIES CONTINUITY EQUATIONS 


The fundamental physical basis for all one-dimensional atmospheric 
models is the set of species continuity equations 


9n . 
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( 1 ) 


where t is the time (sec) and z the altitude (cm) , and is the species 

concentration (cm“^) (molecule units are suppressed). The terms P and L are 
photochemical production (cm“^ sec”^) and loss (sec“^) rates, respectively, 
and (f) is the vertical particle flux (cm“^ sec“^). According to Colegrove 
et al. (ref. 11), we can write the species flux (ignoring thermal diffusion) 
as 
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In (2), K is the eddy diffusion coefficient (cm^ sec~^), and D. is the 
molecular diffusion coefficient (cm^ sec"^) for a species in air. The 
atmospheric scaling factor, y, is a function of the air number density 
M (cm“^). 


M 


(3) 


where is the number density at a reference height, Zq. The scaling 
factor, Y> can also be related to the pressure scale height, E (cm), and 
air temperature, T (K) , using the hydrostatic equation and the ideal gas 
law: 


T 



( 4 ) 


The scale height is a function of temperature. 


E = 



(5) 


where k-^ is Boltzmann’s constant {g cm^ sec”^ , g is the gravitational 

acceleration (cm sec“^) in the lower atmosphere, assumed to be constant with 
height, and m is the average mass of an air molecule (g) ^ equivalent to 
about 29 amu in well-mixed air. As a check on the consistency of our 
adopted atmospheric structure parameters, M and Ty we also calculate y using 
(4) and require the results to be within about 10 percent of those from (3). 
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Moreover, the altitude of the change in the temperature lapse rate in the 
lower stratosphere is correlated with the tropopause height implied by the 
minimum in the eddy diffusion profile. 


The species scale height, (cm) in (2), like is defined by 


H. = 


k 

Bm.g 


( 6 ) 


where is the species mass (g)» Molecular diffusion is only important 

above the "turbopause" at about 100 km. In our model, we use the lower 
thermospheric region between 90 and 120 km as a buffer zone to insulate the 
numerical solutions from upper boundary effects. Accordingly, we have 
adopted one molecular diffusion coefficient for all of the species, and as a 
convenience, we define it by 


D = (7) 

where v-q is a fixed diffusion velocity. Equation (7) implies a dependence 
of D on temperature and pressure, p, which is D ^ T'^lpi this dependence 
is actually very close to reality (ref. 11). We have normalized D to a 
value of 1.5 X 10^ (cm^ sec"^) at 120 km, a typical value for atomic oxygen 
at this height (e.g., see ref. 11). Molecular diffusion is completely 
negligible in the stratosphere — indeed, most stratospheric models ignore it 
altogether. 

When we first began solving atmospheric continuity equations we found 
that by scaling the species concentrations into mixing ratios using y, we 
could effectively reduce diffusion gradients and obtain greater numerical 
stability and accuracy. Hence, in our analysis we use species mixing 
densities defined as 


p . - yn. (8) 

where y can be normalized to any value (we set y = 1 at our lower 
boundary). The species flux equations then become: 


with 
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Accordingly, we can recast the species continuity equation (1) into the form. 
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where and are evaluated using species concentrations, not mixing 
ratios. 

At the outset of our stratospheric studies almost 6 years ago, we 
noted that the solutions of the species continuity equations using approxi- 
mate, but very stable and efficient, numerical integration schemes could 
lead to the violation of mass conservation within groups of related species 
such as the nitrogen oxides (NO, NO 2 , ^ 2 ^ 5 ’ HNO 3 ). We overcame 

this problem by developing a general technique that accurately monitors the 
total concentrations of nearly conserved sets, or families, of compounds 
(ref. 1 ); actually, the concept of aeronomic families has a long history in 
the published literature. Our technique is based on an equivalent continuity 
equation for an entire family of gases which we obtain by summing the 
continuity equations of all the family members after each is multiplied by an 
appropriate factor, a^, representing the weight of the species — or the number 
of odd-atoms of the type considered — within the family. Thus, summing (1) 
for a particular set of species we obtain: 
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where £ is the family index, 


■ ,S, 


(13) 
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and I indicates the subset of species in the Jlth family. 

The family flux, which follows from (9) , is 

a + D) .. ^ . 

y 8a ■ 

The factor r>p in (16) is a weighted average value, 


(15) 


(16) 




(17) 


4 



which is nearly time invariant since the ratios of species concentrations 
within a family are nearly constant under most conditions. Obviously, the 
species boundary conditions, which we discuss later, are also additive for 
families. 

What makes our family technique useful is that the net photochemistry 
for the family, represented by in ( 12 ), usually contains only 

slow processes which affect the total family abundance and precludes the 
rapid chemical interactions that occur between the family members. Accord- 
ingly, the summed continuity equations in ( 12 ) are inherently more stable 
than the individual species continuity equations, and their solution can be 
made very accurate. Another useful aspect of (12) is that these equations 
have the same form as ( 11 ) and, as we will see later, they can be solved with 
the same numerical algorithm. 

Family concentrations obtained from (12) are used to correct the species 
abundances computed using (11). The families of species that we have utilized 
in our model are given in table 1 in the hierarchical order that they are 


TABLE 1.- FAMILIES OF SPECIES 


Hydrogen: 

Hd),"^ 0H(1), H02(1), H202(2), HN02(1), HNOgd) 


Carbon: 

CH(1), CH2(1), CHgCD, CHOd), CH20(1), CHgOd), 
CH 302 d), CHi^ 02 ( 1 ) 


Nitrogen: 

N(l), NOd), N02(1), NOgd), HN02(1), HNOgd), N 2 O 5 
C 10 N 02 ( 1 )^ 

( 2 ), 

Sulfur : 

Sd), sod), S02d), SOgd), HSOgd), H2S0^^(1) 


Chlorine-oxide : 

Old), Clod), 0102 ( 1 ), C10N02(1) 


Odd-chlorine: 

Cl(l), ClO(l), C102d), C10N02(1), HCld) 


Bromine : 

Br(l), BrO(l), HBr(l) 


Oxygen: 

0(1), 03 d), O^Dd), N02(1), N0g(2), HNOgd), 

(3), 


ClO(l), 0102 ( 2 ), C 10 N 02 ( 2 ), BrO(l) 



^The number in parentheses is the species weight within the family. 


Chlorine nitrate is not corrected in the nitrogen family, but rather 
in the chlorine-oxide family. 

solved and applied to correct the species concentrations. For example, HNO^ 
is a member of the hydrogen and nitrogen families, but it is adjusted during 
the hydrogen correction cycle and remains fixed during the nitrogen cycle. 
Two types of families are easily recognizable: one in which the family 
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members recombine in pairs (hydrogen, nitrogen, oxygen), the other in which 
the members recombine with species outside of the family, or do not recombine 
at all (the other families). Families can also be classified according to 
their internal chemistry: the members of some families cycle rapidly among 

themselves in reaction loops (hydrogen, nitrogen, oxygen, chlorine), while 
the members of other families react in a chain, one species leading to the 
next, from the family source to its sink (carbon, sulfur). It is noteworthy 
that the consideration of families forces one to identify the atom carriers 
in aeronomic processes. For example, in the oxygen family we can readily 
find reactions that might not seem to involve oxygen atom transfer or 
recombination, but actually do; thus, the reaction of CIO with NO to form Cl 
and NO 2 does not Involve any odd-oxygen production or loss according to our 
criteria for odd-oxygen (see table 1) . 

The family photochemical terms in (14) can be written in two ways, 
depending on the type of family being considered: 


'SI 


- 


for families whose members recombine in pairs, and 


(18) 


SI 




(19) 


for the other families. These relations are logical extensions of the 
fact that the concentrations of family members are usually fixed fractions 
of the total family concentration, at least for short time periods. 


There are two simple ways of correcting species concentrations using 
the family concentration. In the first, only the most abundant member of the 
family is adjusted to bring into agreement the quantities 


I " ^ ^ 

H}. 


( 20 ) 


In the second method, all of the family members are scaled by the ratio 


W ( 21 ) 

In the following section we discuss our application of the family correction 
scheme to our finite difference solutions for the species concentrations. 


THE NUMERICAL SOLUTION OF THE SPECIES CONTINUITY EQUATIONS 


Before presenting our detailed finite difference analysis of the species 
continuity equations, it is appropriate to describe the basic mechanics of our 
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computer model. The model extends from 10 to 120 km over a 56-point altitude 
grid with a 2-km vertical spacing. We presently compute the concentrations 
of 47 species using 123 reactions and 31 photoprocesses. Not all of the 
species are currently **active" ones in that their concentrations have been 
effectively set to zero for calculations (e.g., ammonia and the bromine com- 
pounds). Data that we initialize in the model are the atmospheric number 
density and temperature, vertical eddy diffusion coefficient, and species 
concentrations. The code internally calculates chemical rate constants and 
photodissociation rates using tabulated photochemical data (e.g., reaction 
rate activation energies and molecular absorption cross sections). 

The time step control is managed as follows: for nondiurnal calculations, 

the initial time step t (which is usually the minimum time step allowed) is 
10^ sec. After each computation cycle, the maximum fractional change in any 
of several critical species (0, 0^, NO, NO^, H, OH, HO^, CH 3 O 2 , CH 4 O 2 , SO^, 

CIO, HCl, BrO) at any height is compared to a preset limiting value e 
(usually 0 . 1 ), and if the change is more than e, t is halved; if the change 
is more than twice e, the computation cycle is repeated until either the 
fractional change is less than 2e or the minimum time step is reached. When 
the fractional change is less than e/2, x is increased by 25 percent. Once 
T exceeds 10^ sec, it is increased by 25 percent only when the fractional 
change per step is less than (e/2) (10^/x) , up to a maximum time step of 
10 ^ sec. For diurnal cycle calculations the time steps are fixed in a 
temporal grid over the day with the number of steps specified as an input 
(usually 100 steps.). The diurnal time increments are calculated by dividing 
24 hr into the required number of steps using a weighting function 

fix) = [ 0 .O 8 + 0.92(|cos xl)^/^] 

where x is the solar zenith angle corresponding to a given time of 
day. 

We write the continuity equations (11) in the finite difference form: 


( 22 ) 

the 
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ik 
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I (23) 


where subscript i is the species index, k is the altitude level, and 
superscript j is the discrete time index. In our notation, quantities that 
are evaluated at the beginning of a time step, j, are "explicit" and are 
known, while quantities which are evaluated at the end of the step, that is. 
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at j + 1, are "implicit" and are to be determined durii^ the course of the 
solution. The altitude increment is 7z(2xlO^ cm), and P replaces yF 
in (11). The diffusive flux divergence is calculated using fluxes centered 
halfway between altitude levels; this assures exact mass conservation for 
vertical diffusion. The terms and y computed at the midlevel 

points, using a logarithmic interpolation, or equivalently, as 


for example. 

Equation (23) can be put in the convenient form, 


'i,k+rik " ^ik "ik " ^i,k-rik "ik 


i = 1, 2, . . ., J (25) 


with 
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ik T i-k 
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7+ 1 

The solution matrix for the in (25) is tridiagonal, and there is a 

simple and fast technique for its inversion (e.g. , ref. 12). The method 
is based on a coupling equation between the solution values at adjacent 
altitude levels, 


P 


J+1 

ik 



- W. 0 ^'^^ 

^ik^ijk+l 


(27) 


Substitution of (27) into (25) leads to the recurrence relationships for 
U and Vi 


ik ■ {ik - i.k-x^kil^ik - i.k-i'^ii 
ik - -hkj{ik - i.k-^'^^i) 


(28) 
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To proceed with our finite difference solution, we must first discuss the 
species boundary conditions in our model. 

For each species, we specify a lower boundary flux and an 

upper boundary flux, which may be fixed or have a specified time 

dependence. But at the lower boundary we also include a flux component. 




( 29 ) 


by defining the ’Velocity” at the boundary, and concentration, 

(if is zero, our code sets <j>^^ = 0 automatically). Whenever we use 

boundary condition (29), we usually set the boundary velocity, V^, to 
1 cm sec"^, which is a typical value for the troposphere. The flux condition 
(29) is equivalent to a solution of the steady-state species continuity 
equation between the ground and the ’ lower boundary at 10 km using known sur- 
face boundary conditions and tropospheric process rates for a particular 
constituent. In other words, equation (29) establishes a loose connection 
between tropospheric processes and the species boundary conditions at 10 km. 
Hence, boundary condition (29) can be used to summarize, in a simplified way, 
the effects on species concentrations of the chemistry and motions in the 
lower atmosphere. Considering the crudeness of our knowledge about tropo- 
spheric aeronomy, a very simple treatment of tropospheric processes in this 
manner is certainly appropriate. 


The use of boundary condition (29) for long-lived, well-mixed gases 
like CO 2 , N 2 O, and CH^ is straightforward —we set to the concentrations 

we wish them to attain, and they are adjusted automatically. The ambient 
troposphere is normally an efficient sink for many gases (e.g., the nitrogen 
oxides, HCl) . In these cases, we usually set wq to a small concentration, 
often to 1, unless observations in the upper troposphere are available. We 
treat the nitrogen oxides somewhat differently, however. First, we specify 
their total lower boundary concentration, 5^(N0^), which we use to calculate 
their total boundary flux; then we apportion this flux among the individual 
species according to their instantaneous abundances. For chemically active 
radicals (e.g., 0, OH, Cl) we simply set ~ ~ also, n-Q = 0. In 

the past we have only exploited the more general application of boundary 
condition (29) during specific pollution studies involving fluorocarbons and 
HCl. 


In our model, when we utilize an eddy diffusion profile with a sharp, 
very stable tropopause level, usually located between 13 and 16 km, the 
boundary conditions at 10 km have little effect on our solutions (except for 
the uniformly mixed constituents, of course). In this situation, the trans- 
port barrier at the tropopause effectively decouples the troposphere and the 
stratosphere for many of the air constituents. 

The boundary conditions that are currently used in our model are 
summarized in table 2. During computer runs, we print out the boundary 
specifications and the boundary fluxes for each species and use them to 
check the accuracy and convergence of our solutions. For example, we can 
balance the hydrogen atom flow via H 2 , H 2 O, and CH^ at 10 km, and we can 
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TABLE 2.- SPECIES BOUNDARY CONDITIONS 


Species 


Lower boundary concentration Fixed upper boundary flux 


at 10 km, cm 


at 120 km, cm 


— 2 —1 
. ^ sec ‘■ 


3.0(11) 


2 . 6 ( 12 ) 


4.0(12) 

3.0(13) 


- 1 . 0 ( 11 ) 


- 2 . 0 ( 8 ) 


1.0(7) 


CF,C1, 

Z ^ 


CH3CI 


1.0(13) 

3.0(11) 

3.0(15) 

1 . 0 ( 6 ) 

1.6(9) 

9.0(8) 

1.0(9) 

5.0(9) 


All species have zero fixed lower boundary fluxes. 

Species not listed here have zero lower boundary concentrations, n^, 
and zero fixed upper fluxes, (})^, except for , CH^ 02 , SO^ , H 2 S 0 ^, and 

CIONO 2 which have = 1 . 

Q 

Inactive species not listed here are: NH^ , Br, BrO, HBr. 

^3.0(11) = 3.0x10^1. 

Q 

The total lower boundary concentration^ for the nitrogen oxides is 

1 . 0 ( 6 ). 

match the total chlorine atom efflux from the stratosphere against the 
integrated production rate. 

With the boundary conditions so defined, we can not proceed with the 
development of our numerical solution. The species continuity equations at 
the lower boundary are written: 





( 30 ) 



- - 




where we recognize that the lower boundary- layer thickness is only h/2 cm. 
The flux condition at the lower boundary is expressed as 


1.0 




- 


.j+i _ j+i . 


- 


1)P 


J+1 

il 


(31) 


In our scheme we have allowed for an ” image” concentration, p . , a distance 

7 + 1 

h below the lower boundary. After eliminating p. from (30) and (31), the 
resulting expression can be rearranged to yield 


^2 


(32) 


Equation (32) defines and for each species, and these quantities 

can be used with (26) and (28) to determine the U and V terms at every 
other altitude level for all the species. 


We note that our formulation of the boundary equations allows the species 
boundary concentrations to be affected not only by local flux divergences, 
but by photochemistry as well. A simpler, well-known boundary specification 
which requires that 


(|) (boundary) = 






- V 


D 


(r. - l)p 


J+i 

il 


(33) 


artificially connects the solution for p.^ and together by a mixing 

ratio gradient attributable only to transport. 


The upper boundary condition is treated like the lower one. That is, we 
write a continuity equation at the upper level analogous to (30) except that 
we utilize an ”image” density above the upper boundary, , (where u is 

the index for the uppermost altitude interval). We also have a flux condition 
like (31) for the upper boundary. Moreover, we know the recurrence relation- 
ship, (27), fqr level k=u- 1 (i.e., we have previously determined the 
quantities and ). Accordingly, with these three equations 

1 . 3 W— 1 “Z-j Zf— J- 
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we can eliminate 
be: 


p% ^ and p . and solve for p . , which turns out to 


p 


J+1 

iu 


bP. - 2(y A)<#). - 2cP. U^. , 

VU U ^U 1 


(1/T) + 4 f - 1) - 2rf^(i + 


(34) 


The species concentrations below the upper boundary are obtained at 
the descending levels, u - 1, u - 2, ... 1, using (27), the known V and 7 
coefficients, and the calculated species concentration at the next highest 
level. Obviously, by solving the species continuity equations in sequence, 
i = 1, 2, . . . , J, only one pair of U and 7 vectors are needed, and these 
can be recycled for each species. 


The family continuity equations described by relationships (12)~(19), 
are solved in the same way as the species continuity equations. There are 
two procedural differences that must be mentioned, however. First, the 
family photochemical loss rate in (12) is linearized in the implicit variable 

which leads to the following numerical generalization of (18) and (19): 



Ik Slk 


(35) 


The definitions of Pq and which are straightforward, are given by 

Turco and Whitten (ref. 1). The second procedural difference is that the 
family boundary conditions must be established by summing (30) and (31) — 
and the corresponding relations at the upper boundary — for each family. 
During this exercise, all the terms that include a species subscript 
(i.e., those containing ^i:u) redefined as 

average values at the beginning of the time step; for example. 



= E 


7 

a .V .p 





(36) 


The family continuity equations (12) and (16), when cast in the finite 
difference form of (23), and after the application of the specialized treat- 
ment of the photochemical terms and boundary conditions just described, 
are solved algebraically exactly like equations (23)- (34) for the individual 
species . 

For our model, the extra computer time required to handle the families 
is about 10 percent of the total computer expenditure, which is more than 
offset by a large gain in numerical speed and stability. Another numerical 
note: since we determine the species concentrations first, we can then use 

7 7 i 

both the , or their average value, to compute the family 

chemical production and loss rates, thereby further improving the solution 
characteristics . 
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The family concentrations are used to correct the computed species 
concentrations in our model. In this sense, our method is more sophisticated 
than techniques that only compute family concentrations and then assign each 
member an abundance based on steady-state photochemical ratios. Our basic 
correction scheme, which is performed automatically in the code, is to adjust 
only the most abundant family member; this leads to the correction equation 


P 


J+1 

ik 


(corr. ) 




(37) 


For the carbon and sulfur families, however, all of the species are 
adjusted using a scaling factor. 




ik 



(38) 


Our correction scheme is effective for at least two reasons: 

1. In most cases, numerically generated concentration imbalances 
among the members of a family are quickly redistributed by the fast 
photochemical coupling among these species. 

2. During a calculation, the species concentrations are restricted 
to maximum changes of about 10 percent or less per step, and the 
family concentration variations are usually much smaller than this, 
so that the actual correction per step is normally a very small 
fraction. 

We note that a species is adjusted only once in the correction hierarchy 
given in table 1. In some cases, at specific altitudes, the correction 
scheme itself may become unstable because the numerical system is over- 
determined. For example, when two species in a family have almost equal 
concentrations, an unstable oscillation can develop where one species is 
corrected at one time step, the other at the next step. This problem is 
eliminated by explicitly stating which species is to be corrected each time. 
These isolated instabilities are easily handled in the computer program. In 
our code, we have achieved absolute solution convergence with time steps as 
large as 10^ sec (i.e., about 3 steps per year), and we could easily do even 
better. 


AERONOMICAL MODEL PARAMETERS 


We will now discuss some of the aspects of the atmospheric physics 
and chemistry which we have included in our model. 
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The rate coefficients used in our calculations are temperature and pres- 
sure dependent, wherever appropriate. For the most part, we have adopted the 
standard rate constant values from published data lists (e.g. , the National 
Bureau of Standards Reports — also see refs. 4 to 10 for tabulations of our 
rate coefficients). Whenever a new rate constant is introduced in our work 
or an old one is adjusted to achieve some aeronomical state, it is usually 
discussed in detail. For several important altitude-dependent rate coeffi- 
cients, most notably those for the formation reactions of HNO^ from OH and NO^ 
and of N 2 O 5 from NO 2 and NO 3 , we have developed analytical pressure- and 
temperature-dependent expressions for the coefficients using steady-state multi- 
level kinetics, thereby eliminating the need for tables of data (see ref. 9 ). 

The rate constants for the critical reactions, 

OH + HO 2 H^O + O 2 (39) 

0 + HO^ OH + 0^ (40) 

are still controversial. In our model, using values of 5.0x10“^^ cm^ sec”^ 
for and 3 . 0 x 10 ”^^ cm^ sec“^ for leads to a total integrated 

daytime OH column abundance of about 8x10^’^ cm , close to that detected 
by Burnett (ref. 13). It is noteworthy, however, that the variability in 
Burnett’s data encompasses nearly an order of magnitude about a value of 
1x10^^ cm"^. Our stratospheric OH concentrations are still somewhat lower 
than the observations of Anderson (ref. 14) after adjusting his values 
upward by a factor of 2 to roughly account for the ratio of OH abundances in 
full daylight to those at the experimental zenith angle of 80 °. 

In our computer code, rate constants and the corresponding chemical 
production and loss terms in the species continuity equations are manipulated 
automatically; a five-number label for each reaction is used to specify the 
reactants and products. Moreover, by scanning the reaction label and com- 
paring the species involved with templates of family constituents, the 
reactions affecting the families are flagged for quick future reference. As 
a check on the automatic chemistry routine, the atom balances for each 
reaction are displayed as output at the start of every computer run. 

We calculate atmospheric photodissociation rates using the well-known 
optical depth formulations of light absorption, and the set of wavelength 
intervals listed in table 3. In this case, 

J = ^ F o exp[-T (a,x)] (41) 

V ^ IJ VP ^ P 

where c7^ is the photorate (sec“^) for process v; p is the wavelength 
interval index; is the total incident solar flux (photons cm~^ sec”^) 

at the top of the atmosphere in the wavelength interval p is the cross 

section (cm^) for process v in interval p; and is the corresponding 

optical depth at the height z and solar zenith angle x* The optical depth 
is defined by 
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TABLE 3.- WAVELENGTH INTERVALS USED TO CALCULATE PHOTODISSOCIATION RATES 


Wavelength range, 
nm 

Number of bins 
in the range 

Bin width, 
nm 

121.6 

1 

0.1 

137.5-177.5 

8 

5.0 

177.5-177.75 

1 

.25 

177.75-201.25^ 


.5 

201.25-202.5 

1 

1.25 

202.5-347.5 

29 

5.0 

347.5-350 

1 

2.5 

350-750 

8 

. 

50 


^In this wavelength region, each 0.5-nm bin is subdivided into 
five 0.1-nm intervals for 0^ Schumann-Runge band absorption data. 


SR 


T^(3,X) = C2(3,x)a2^ + 


(s,x) + 


(42) 


where subscripts 2 and 3 refer to O 2 Q^nd O 3 , respectively; C in each 
case is the integrated molecular column (cm from the point of observation 
(s,x) to the sun, and are continuum absorption coefficients (cm^) , and 

is the 0^ Schumann- Runge band optical depth, which we will discuss shortly. 
We have ignored other atmospheric absorbers of solar radiation such as NO^ 
since their effects on photorates are negligible. 


Our solar fluxes are taken from Ackerman (ref. 15) above 300 nm and 
Donnelly and Pope (ref. 16) below 300 nm; the latter fluxes, for a moderate 
level of solar activity, represent a compromise between Ackerman’s larger 
values and recently observed lower values (e.g., ref. 17). In our calcula- 
tions we do not account for Rayleigh multiple scattering of sunlight, but 
we do roughly account for the effective planetary albedo by increasing the 
incident flux, F, by an albedo factor, a (e.g., see ref. 18): 


so that 


a 



X > 320 nm 
300 < X < 320 nm 
X < 300 nm 


F ^ ^(1 + a) 


(43) 


(44) 
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We calculate the integrated O 2 column density, using the optical 

depth factor approximation for an exponential atmosphere which is based 
on the Chapman function (refs. 19 and 20). For this calculation we assume 
that O 2 comprises 21 percent of the total number of air molecules, and we use 
our model temperature profile to evaluate the appropriate scale heights. We 
calculate the ozone column, C 3 , by numerically integrating the instantaneous 
ozone distribution in the model along a ray from the point of observation to 
the sun using a logarithmic interpolation scheme (ozone concentrations below 
10 km are specified as input, and those above 120 km are extrapolated 
exponentially assuming a 3-km scale height). Thus, 

C3 = ^ 

k 




fOo], 


(45) 


where is a path length defined by 


s 


k 




(46) 


with being the distance (cm) from the Earth’s center to the altitude 

level k, and the perpendicular distance (cm) to the (extended) ray of 

observation, and where 


with the brackets ([ ]) indicating a concentration. 

Before Hudson and Mahle (ref. 21) had published their parameterized 
equations for calculating 0^ Schumann-Runge band absorption, we had already 
treated this absorption in detail using a band model. Blake et al. (ref. 22), 
who collected 0^ absorption data in the S-R band region with an instrument 
having about a 0.1-nm resolution, showed that the observed 0^ S-R band 
optical depth depended simply on the square root of the 0 ^ column density for 
a wide range of optical depths. This square root absorption law has been 
used by Brinkmann (ref. 23) to study water vapor dissociation in the terres- 
trial atmosphere. In our model, we define the 0^ optical depth in the 
Schumann-Runge bands as 
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( 48 ) 


Here is the 0^ absorption cross section obtained from Blake et al. 
(ref. 22) at 0.1-nm intervals in the S-R band system (177.75-201.25 nm) . 
In analyzing the Blake et al. data, we have accounted for 0^ continuum 
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absorption and its apparent increase at high experimental absorption cell 
pressures (see ref. 2). The parameters and are the weak and 

strong absorption limits beyond which the simple square root absorption law 
is inapplicable. In the weak absorption limit, we choose = 0.15 to 
obtain the correct integrated oscillator strength for the bands. Blake et al. 
(ref. 22) found that their absorption cross sections displayed a square root 
dependence on the column density up to ^ 3. In numerical experiments, 

we have found that none of our photorates are influenced by more than 5 percent 
for any selection of > 2; accordingly, we let ^ “ in our model. 


Using a band absorption model, the 0^ photodissociation rate in the 
S-R bands is easily calculated since it is simply proportional to the 
derivative with respect to of the corresponding optical transmission 
factor. In (41), this is equivalent to substituting for 


SR 

o 



(49) 


where 



SR ^ 

T < 


SR 

w 


or 


SR SR 

T > T 

S 


SR 

w 


< T 


SR 


< T 


SR 

s 


(50) 


Except for the 0^ S-R band data, the absorption cross sections for other 
species are specified every 0.5 nm in the S-R band region (although they 
probably only need to be specified every 2.0-5. 0 nm) . In each of these 
0.5-nm wavelength bins, we calculate an average 0^ S-R band transmission 
factor and dissociation cross section: 



e 



(51) 



(52) 


where the index g ranges over the five 0.1-nra bins within the 0.5-nm 
interval. The quantities defined by (51) and (52) are inserted directly 
into (41) . 

We have made detailed comparisons (refs. 2 and 24) between the O 2 
absorption and dissociation profiles calculated with our band model and those 
from the Hudson and Mahle model (ref. 21); the agreement is better than about 
10 percent in the stratosphere and about 25 percent in the mesosphere. 

This correspondence is excellent, considering that the two results are derived 
from completely independent data bases, and that no tuning of the band model 
(to account for temperature effects and instrumental uncertainties on the 
observed low-resolution absorption cross sections) has been performed. The 
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band absorption model also preserves the simple optical depth formalism and 
is therefore easy to include in existing photodissociation rate algorithms. 

Our species photodissociation cross sections are given in reference 2. 
Newer data, not listed in reference 2, are used for the 0(^Z}) quantum yield 
from 0^ photolysis (refs. 25 and 26) and for the absorption cross sections of 
N^O (ref. 27), HCl (ref. 28), CIO (private communication from Jaffe, Ames 
Research Center, 1976), and CIONO^ (ref. 29). 

In our model we can select one of several options for computing 
photorates. For diurnal calculations, the rates are calculated at a fixed 
set of zenith angles corresponding to the fixed diurnal time grid mentioned 
earlier (the latitude and season which are selected also determine the zenith 
angles). For nondiurnal runs, we can compute photorates at a fixed zenith 
angle specified by a latitude, earth declination, and time of day. Or, we 
can compute, in two different ways, photodissociation rates which are 
averaged over a 24-hr day: 1) by determining the rates at up to 50 zenith 

angles over one-half the day and numerically averaging these values; and 
2) by estimating the rates using an approximation described by Rundel^ for 
spring/fall conditions: 


J - 


c 

o 


1 + b 

o 


J(x 


noon 


) 


(53) 


with = 0.5, = 0.835, = 0.619, and Xnoon ^ have found 

that Rundel’s approximation gives results in good accord with the diurnally 
time-averaged rates. 

Cogley and Borucki (ref. 30) have calculated average photodissociation 
rates using the approximation 


J = (54) 

We have compared photorates computed using (53) with those from (54) when 
5^ = 0.5 and Xq ^ 55°. For the most part, the values are within 5-10 per- 
cent of each other. However, in regions of strong absorption the differences 
become quite large, and in this case, Rundel’s predictions are generally 
more satisfactory. Still, it is not clear whether, in the O 2 Schumann-Rung e 
band absorption region, more than one set of parameter values is needed for 
the Rundel model because of the complex dependence of the S-R optical depth 
on the molecular oxygen column density. 

Modeling techniques that only use averaged photodissociation coefficients 
to simulate diurnal effects can often give incorrect and misleading results 


^Private communication from R. D. Rundel, Johnson Space Flight Center, 1976. 
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for the abundances of many stratospheric gas constituents and for the net 
depletion of ozone by chlorine oxides. Therefore, we have developed a 
simple scheme for modifying the species continuity equations that accurately 
accounts for the effects of diurnal variations on computed species concentra- 
tions and ozone perturbations (ref. 3). To do this, we average the continuity 
equations over a 24-hr day and derive equivalent relations for the average 
concentrations n • : 


. 

% 


= E 

Qi'k- 


e . 3 -7 .7 -7?- 


d 

3a 


(f) (n . ) 


(55) 


Here, is a generalized rate constant which can be a photodissociation 

rate; is ±1 depending on whether the species i is produced or 

destroyed, respectively, by reaction (j/c) ; and is the average species 

flux evaluated using the height distribution n^. In our approximation, the 
average species concentration is divided into two components, the 

average daytime and nighttime species abundances rP. and rP. ^ respectively, 
so that : ^ ^ 


T 

D "-D ^ 
n . = n . ~ + 
^ ^ T 


T 

N ^ 
i T 


n 


(56) 


where T-q and Tj^ are the daytime and nighttime durations, respectively, 
and T is the sum of T]j and During a given diurnal cycle, a species' 

daytime and nighttime concentrations are assumed to be constant at their 
respective average values, rp and rp \ hence, our solution scheme parameterizes 
diurnal variations in the form of a simple, two-level step function. With 
this assumption, it is then convenient to define the quantities 


a 




T 

N 

T . 

^ T 


(57) 


where is the nighttime-to-daytime species concentration ratio 


r . = 


rP. /n 

7 ■ 


(58) 


In terms of the parameters 
6., in (57) is 



and a^. 



the rate constant scaling factor 



(59) 


The 3 factor is obviously 1 unless both of the reactants have diurnal 
variations. In (55) the photodissociation processes take the form 


.n ./a . 

V y % 


( 60 ) 


19 



where J is a 24-hr averaged photorate that can be computed using any one 
of the several techniques discussed earlier. 

Our diurnal averaging scheme quantifies the effects of nighttime 
chemistry on the average day-to-day abundances of air constituents by adjust- 
ing the appropriate photochemical rate constants to account for the presence 
at night of certain reactant gases. Our step function approximation for 
diurnal variations greatly simplifies the evaluation of the 24-hr average 
rates of the nonlinear photochemical interaction terms in the species con- 
tinuity equations. The relative magnitudes of the two-level diurnal species 
abundances are parameterized using the night-to-day concentration ratios. 

For most stratospheric constituents, these concentration ratios are 
either effectively 0 or 1. For other species we can calculate r values by 
using simplified nighttime chemical reaction schemes, assuming as initial sun- 
set conditions the appropriate daytime species concentrations and ignoring 
the effects of transport on the ratios. This approach leads to t values 
that are functions of the model (rate constants and species concentrations) , 
and can be updated accordingly during computer simulations. Actually, the 
concentration ratios are nearly invariant quantities that could, for example, 
be calculated using the results of a diurnally varying model simulation. 

The averaged species continuity equations in (55) have exactly the 
same form as (1) and can therefore be solved in the manner we have already 
described. In fact, in our model we can turn our averaging procedure on or 
off quite simply. Furthermore, since the form of the family continuity 
equations is also unaffected by our averaging technique, their solution also 
proceeds just as before. Once the average species concentrations, n^, are 
calculated, we can obtain daytime and nighttime abundances using (57) and (58). 


We have compared the predictions of our averaged model with those of a 
diurnally varying model (ref. 3); typically, the calculated concentrations 
are within 10 percent of each other, but with larger differences of up to 
20 percent occurring at some altitudes (usually where a species is decreasing 
rapidly and has a very small absolute abundance). This is in striking con- 
trast to the large errors obtained from a model using only averaged photodis- 
sociation rates where, for example, an order of magnitude underestimate of 

the N^0_ concentration results. 

Z o 


Our model calculates the light emission intensities of vibrationally 
excited OH^ and electronically excited 02 (^Ag, ^^g) molecules in the 
atmosphere. Predicted column- integrated emission rates have been compared 
with observations to provide a constraint on the model parameters. For OH^ , 
photon efficiencies (i.e., the number of photons emitted for each OH^ formed 
by reaction) in the vibrational sequences At? = 1, 2, 3 are calculated with 
a detailed vibrational state model for 0H(Z ^11, t? = J , 9), which includes 
radiation cascading and chemical and collisional quenching. For the 0^ 
singlet delta states, we have incorporated all of the known photochemical 
excitation and quenching mechanisms in our computer simulation. 
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We have already mentioned our eddy diffusion coefficient, Ky whose 
values are given in figure 1- Our profile is similar in its general 
characteristics to one originally 
proposed by Wofsy and McElroy 
(ref. 31), but with some important 
differences- Our dif fusivities 
are based in part on meteorologi- 
cal considerations and in part on 
the successful prediction of the 
observed distributions of several 
atmospheric tracers *— in particu- 
lar, methane, nitrous oxide, and 
carbon 14- The troposphere, 
which is usually marginally 
stable or unstable, is usually 
subjected to strong mixing motions 
that are characterized by a large 
vertical eddy diffusion coeffi- 
cient (-3x10^ cm^ sec“^). 

However, above the tropopause, 
where the temperature lapse rate 
decreases sharply, the strato- 
spheric air mass is highly stable 
and very resistant to vertical displacements. Typically, vertical velocities 
just above the tropopause are only about 1 percent of those found in the 
troposphere. Hence, we expect the eddy diffusivity also to change abruptly 
at the tropopause, as indicated in figure 1. In our model we have fixed 
the tropopause height at 13 km, which is a representative value for midlati- 
tudes. Even though the entire stratosphere is thermodynamically stable, 
we know, by observing the distributions of inert tracers like carbon 14 
in the lower stratosphere and photochemically active species like methane 
at greater altitudes, that the effective vertical diffusion velocity increases 
with increasing height above the tropopause. This property is reflected 
in the shape of our eddy diffusion profile above the tropopause. 

The stratospheric distributions of carbon 14 resulting from atmospheric 
thermonuclear explosions have been measured following the nuclear bomb tests 
of the early 1960s (e.g., see ref. 32). Accordingly, we have adjusted the 
lower portion of our eddy diffusion profile so that our model predicts 
distributions which, even after 2 years of atmospheric relaxation following 
a nuclear detonation, are still close to the observed values. This agreement 
is illustrated in figure 2(a). Johnston et al. (ref. 32), in their analysis 
of the carbon 14 data, have concluded that it would be difficult to reproduce 
most of the carbon 14 observations using eddy diffusion coefficients with a 
tropopause level higher than 13 km, and we also find this to be true. 

In the upper stratosphere, the predicted vertical distributions of 
selected ” tracer” gases such as methane and nitrous oxide are not independent 
of photochemistry. Indeed, new laboratory measurements of photolysis cross 
sections and reaction rate constants often necessitate adjustments of the 
eddy diffusion profile. In addition, the selection of eddy dif fusivities 



DIFFUSION COEFFICIENT, Cm^/seC 

Figure 1.- The eddy diffusion coefficient 
{K) profile used in some of our one- 
dimensional model studies. 
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ALTITUDE, km 



AVERAGE EXCESS AT 30® N IN 10^ 


atoms cm"^ 

(a) Carbon 14 observations are rep- 
resented by solid lines (see 
ref. 32), predicted distributions 
by broken lines. 


at the higher altitudes also depends 
on the manner in which constituents 
with large day-to-night concentration 
variations are averaged over the 
diurnal cycle. For example, methane 
is decomposed principally by its 
reaction with OH in the stratosphere 
and mesosphere, 

OH + CH^^ ^ H^O + CH 3 (61) 

Since OH disappears rapidly after 
sunset, as a consequence of reaction 
(39), the predicted total amount of 
methane destroyed in a day is less 
for a model which diurnally averages 
the OH abundance than for a model 
which does not. In our case, using 
the averaging procedure outlined 
earlier, the upward flux of CH^ 
required to maintain its high-altitude 
distribution is smaller than would 



(b) Observed and predicted vertical 
distributions of methane and 
nitrous oxide. 


Figure 2.- Vertical distributions 
of carbon 14, methane, and 
nitrous oxide. 


otherwise be expected. Therefore, 
our diffusion coefficients have been 
reduced somewhat above 30 km to 
compensate for this effect. Our 
predicted methane and nitrous oxide 
concentrations, shown in figure 2 (b), 
are in good agreement with the lim- 
ited observations that are available. 

We have recently developed an 
aerosol particle model to complement 
our one-dimensional photochemical 
model. The aerosol physics treated 
in the model, and the numerical proce- 
dures required to solve the related 
aerosol continuity equations, are 
described in detail elsewhere (ref. 10 ). 
In the model, SO 2 and OCS molecules 
diffuse upward from the troposphere 
into the stratosphere; sulfur gases are 
also injected there by volcanoes (and 
aircraft). At high altitudes OCS is 
photolytically decomposed by ultra- 
violet radiation, and the products 


quickly react to form SO 2 (ref. 33). 

Sulfur dioxide is subsequently oxidized into sulfuric acid vapor via the 


intermediary radical HSO 3 , which is assumed to have a short chemical lifetime 
in air. Once formed, H 2 S 0 ^ molecules nucleate with water vapor onto 
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condensation nuclei which are also transported upward from the troposphere. 
The resulting acid solution droplets continue to grow by heteromolecular 
condensation of water and acid vapors, coagulate with one another, settle 
gravitationally, and diffuse by turbulent mixing. Droplets that rise above 
the aerosol layer (to about 30 km) evaporate rapidly to their bare cores, 
which settle again to lower altitudes. Evaporation is a source of 
molecules above the aerosol layer, but this does not materially affect the 
shape of the layer. Our model also follows the accumulation and evolution 
of the solid cores within the aerosol droplets. 

Preliminary calculations with the aerosol model are in excellent agree- 
ment with numerous observations of the natural particulate layer. The com- 
puted peak mass mixing ratio for the layer is near 21 km, some 8 km above the 
tropopause, as observed. Particle number mixing ratios are close to measured 
values, and the droplet size distribution and composition are also correctly 
reproduced. 

Several feedback mechanisms are built into our atmospheric model. Of 
course, the photochemistry and parameterized dynamics are completely inter- 
active, which significantly influences the response of the model atmosphere 
to simulated perturbations. Aerosol particles can also interact with the 
gaseous compounds in our simulation, although we have not yet included possi- 
ble radiation and surface catalysis feedback effects. The ozone profile in 
our model controls, to a large extent, the ultraviolet radiation fluxes 
reaching the stratosphere. Accordingly, we recalculate all of the photodis- 
sociation rates (and rescale the rate constants during diurnally averaged 
runs) whenever the integrated ozone column above 10 km, or above 30 km, 
changes by a specific fractional amount. This fraction is usually selected 
to be small relative to the expected ozone variation in order to achieve good 
numerical resolution without having to recalculate the photorates at each 
time step, which is very inefficient. 

We have developed a stratospheric temperature model that utilizes the 
ozone ultraviolet heating rates of Lacis and Hansen (ref. 34) and the CO 2 
infrared cooling rates of Dickinson (ref. 35). The air temperatures are 
recalculated whenever the ozone abundance changes significantly, as described 
above. The atmospheric scale heights and densities and the rate constants 
are also redetermined at these times. The temperature subroutine is pres- 
ently inactive in our model since we discovered by experimentation that 
temperature feedback effects on ozone concentrations were only about 
10 percent or less of the overall ozone variation (resulting in slightly 
smaller predictions of ozone reductions in polluted air). Moreover, we do 
not feel that parameterized treatments of heating and cooling that neglect 
the related effects on atmospheric dynamics are valid for the large ozone 
perturbations often studied with our model. 

Finally, we mention another unique feedback mechanism in our model 
which is related to the use of families of species. Our numerical solution 
of the family continuity equations is largely independent of the solution 
for the individual species distributions. However, in order for the complete 
computation to be stable and to proceed to a steady state, these two distinct 
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solutions must converge to the same species concentrations. Hence, in our 
code we have a built-in detector of numerical inconsistencies. 


OPERATIONAL MODES FOR MODEL CALCULATIONS 


Our computer model is organized so that we only need to flag several 
input parameters in order to control the mode of calculation and the printing 
and plotting of output data. The choices for computing photodissociation 
rates, and the optional diurnal averaging scheme, have been discussed 
earlier. The basic temporal modes of operation are: 

1. Steady-state runs which extend for a specified time duration. 

2. Time-dependent calculations which are similar to steady-state runs 
except that certain parameters may have explicit time variations, and data 
at intermediate times are of interest. 

3. Diurnal simulations where, in order to achieve rapid solution 
convergence, we usually use diurnally averaged steady-state species profiles 
as initial conditions. 

In any of these temporal modes, we can add physical perturbations to the 
system. The most important of these are: 

1. Supersonic transport exhaust: We inject SST effluents (NO, H 2 O, 

SO^) uniformly in hemispherical or global shells, the amount and altitude 
profile depending on the engine, traffic, and flight path models. 

2. Space shuttles: We add HCl and NO from the shuttle launch rockets, 
spread hemispherically or globally, with an altitude profile characteristic 
of these vehicles at a rate determined by the launch schedule. 

3. Fluorocarbons: Starting with a steady-state ambient atmosphere, 

we inject CF^Cl^ and CFCl^ at the ground at globally averaged rates (which may 
remain constant or change with time, and may end abruptly). Each fluorocarbon 
is allowed to accumulate in the troposphere at a rate proportional to the 
difference between its input flux at the surface and its escape flux into 
the stratosphere (or into the lower model boundary at 10 km). Thus, we 
solve the time-dependent growth equations for the tropospheric fluorocarbon 
content, and we couple these solutions to the one-dimensional model through 
the boundary conditions at 10 km. As we have always done, we allow for the 
possibility of tropospheric loss mechanisms for CF 2 CI 2 and CFCI 3 by assigning 
them average tropospheric lifetimes (from 30 years to <») . 
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CONCLUDING REMARKS 


We have presented a complete and concise outline of our one-dimensional 
atmospheric computer simulation. Some of the unique aspects of our model 
which we have described above are: 

1. A technique for obtaining rapid and accurate solutions of species 
continuity equations using the concept of conserved families of aeronomically 
related compounds. 

2. An averaging scheme for simulating effects of diurnal variations 
on atmospheric constituent concentrations. 

3. An efficient O 2 Schumann-Runge band absorption model for computing 
molecular photodissociation rates. 

4. A detailed simulation of aerosol particle formation and evolution 
in the Earth’s stratosphere. 

As we mentioned in the Introduction, specific examples of our model 
predictions can be found in the literature. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif. 94035, April 14, 1977 
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